function ve = solve_eq(D, c,s,param,glob,options)

    param.D = D;
    v =  solve_valfunc(c,s,param,glob,options);

   v1     = reshape(v.v1, glob.n(1), glob.n(2));
   v1     = griddedInterpolant(glob.B, glob.S,v1, 'linear');   
   ve     = v1(repmat(param.b0, size(glob.agrid)), glob.agrid);
   ve     = sum(glob.a_stat.*ve);

    
end